Scalar field confinement as a model for accreting systems 
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We investigate the possibility to localize scalar field configurations as a model for black hole 
accretion. We analyze and resolve difficulties encountered when localizing scalar fields in General 
Relativity. We illustrate this ability with a simple spherically symmetric model which can be used to 
study features of accreting shells around a black hole. This is accomplished by prescribing a scalar 
field with a coordinate dependent potential. Numerical solutions to the Einstein-Klein-Gordon 
equations are shown, where a scalar filed is indeed confined within a region surrounding a black 
hole. The resulting spacetime can be described in terms of simple harmonic time dependence. 



I. INTRODUCTION 

Self-gravitating scalar field configurations have been 
very useful in many aspects of gravitational theory. Their 
role as describing matter models (eg.P, [3, d, 0, [l]); 
as governing mechanisms to model inflationary scenar- 
ios (eg. [1, 0|); as probes of strong curvature regions 
(eg.[3, Q), etc. has made them an ideal tool in a number 
of fronts. In this work we examine exploiting scalar fields 
to mimic some salient properties of accreting black hole 
systems. To this end, it is desirable to explore a configu- 
ration where the scalar field simulates an accretion disk 
surrounding a black hole. For this purpose, one should 
be able to confine the scalar field within some compact 
region surrounding the black hole. Since massless scalar 
fields radiate away to infinity, the model sought after 
should include a mechanism that will prevent this from 
happening, at least to some non-trivial extent. (The ex- 
istence of bound states for particular cases in spherical 
symmetry are studied in [lO|, ) . 

One way to confine the scalar field would be by em- 
ploying a potential well which would introduce some 
sort of barrier and thus allow for confinement. The use 
of carefully chosen potentials is common practice with 
scalar fields, and are usually functions of the field it- 
self. Examples of this kind of potential are the quadratic 
(F(0) oc (jP) -that introduces a mass term-; and the 
quartic {V{4)) oc (j)^). However, that kind of potential 
does not allow for confining the scalar field within a spe- 
cific region of space, that one can specify a priori. 

What we are looking for is a potential that somehow 
depends on the coordinates and in particular can be cho- 
sen to describe a potential well within a region. How- 
ever, this proposition seems a priori at odds with main- 
taining covariance. The difficulty one encounters with a 
coordinate-dependent potential, is that the correspond- 
ing stress-energy tensor is in general inconsistent, in the 
sense that its divergence will not be zero for a non-trivial 
scalar field. This fact, together with the Einstein equa- 
tions, would imply that the Bianchi identities are not 
satisfied. 

Faced with this situation a possible way of confining 
the scalar field would be to introduce a background with 



respect to which coordinates could be defined. This ap- 
proach would be in line with bi-metric theories of gravity 
(eg-[l3)- Another approach would be to fix a suitable co- 
ordinates already at the level of the action as is done in 
through some suitably introduced Lagrande multipli- 
ers. This procedure then provides a way to convariantly 
adopt coordinates which could, in principle, be used in 
the potential. However, this method would strongly link 
the adopted gauge with the type of potential introduced 
and it is not yet clear whether it can be made of practi- 
cal use. An alternative way, which is the one we pursue 
here, is to exploit symmetry considerations without re- 
sorting to introducing any other feature in the problem. 
The existence of the symmetry provides a simple way to 
consistently introduce a coordinate-dependent potential 
in the problem. Certainly, while more restricted than 
other possible viable options, this approach is the more 
direct one. A particular case of a coordinate dependent 
potential has already been implemented in [3 EBl to 
effectively simulate angular momentum in spherical or 
axial symmetry. 

In this work we concentrate mainly on the case of 
spherical symmetry, but give prescriptions for the im- 
plementation of potentials in both spherical and axial 
symmetry. We will see that, if the space is spherically 
symmetric, we can implement a potential that depends 
on the areal radius. In the same way, for an axially 
symmetric spacetime, the potential can depend on the 
length of the closed integral curves defined by the asso- 
ciated killing vector. Even though one will not be able 
to specify the potential as an arbitrary function of any 
coordinate, one may still be able to confine a scalar field 
to some region, as is it shown in this work for the case 
of spherical symmetry. This fact will become apparent 
in section III Al and in its applications in the rest of this 
work. 

This paper is organized as follows. In section [TTl we 
study the specification of a stress-energy tensor for a 
scalar field with a coordinate dependent potential. Show- 
ing that such implementation is possible when the space- 
time possesses a symmetry. In particular, the case of 
spherical symmetry is studied in depth (we also consider 
an axi-symmetric case in an appendix). In section HITl we 
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describe the formulation used, and the resulting equa- 
tions. In section [IV] we discuss how the equations are 
solved numerically, after obtaining initial data by two 
different methods. In section |V] we show and analyze the 
numerical solutions obtained, finding that, after some 
transient behavior, the scalar field reaches a state de- 
scribed by a simple harmonic time dependence and re- 
mains confined to a region surrounding the black hole. 
We have confirmed these for initial masses of the scalar 
field up to 50% of that of the black hole. Finally, we 
make some final remarks in section IVTl In all this work 
we use Einstein's index notation and geometrized units. 



II. SCALAR FIELD ON A 
COORDINATE-DEPENDENT POTENTIAL 

In this section we study the specification of a stress- 
energy tensor for a scalar field with a coordinate depen- 
dent potential. Our motivation is to somehow confine 
a scalar field within a region around a black hole. The 
resulting system would share features of a black hole in- 
teracting with an accretion disk. We will see that this can 
be done when the space-time posses a symmetry. How- 
ever, the specification of such potential is not completely 
arbitrary since it must depend on the coordinates only 
through some particular function. Knowing the approxi- 
mate dependence of that function on the coordinates, one 
can then construct a potential that confines the scalar 
field. 

Before presenting our approach, we include an 
overview of how the equations of motion are obtained 
from a stress-energy tensor in the case of a coordinate- 
independent potential. Then, based on that proce- 
dure, we will study the generalization to the case of a 
coordinate-dependent potential. 

The equations of motion for a real scalar field cj) on 
a coordinate-independent potential can be derived from 
the stress-energy tensor 



common factor, 



ra6 = rWa6 + r(^)a6 



(1) 



where, for later convenience, we have split this tensor 
into what we call the "kinetic" and "potential" terms: 
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rC^^ab = (Va0)(Vb0) - -<?a6(Ve(/')(V^0) 



(2) 



(3) 



The kinetic part, T^''\ij, corresponds to a massless scalar 
field without a potential. 

The equations of motion can be obtained P, [l^ 
through the condition 



(4) 



which must be satisfied to be consistent with a covariant 
theory. Equation ^ can be re-expressed with Wb<p as a 



= WaT\ = (Vfc0) £(0), 



(5) 



where £((/>) contains second order derivatives of (j). The 
equations of motion for a non-trivial scalar field is then 



£(0) =0. 



(6) 



For example, for 1^(0) = m'^cjP we obtain the Klein- 
Gordon equation. 



£(0) EE (VaV" - m^) 4> = 



(7) 



This is analogous to the Lagrangian approach, where 
the variation of the action is set to zero, and, after 
integrating by parts, the integrand becomes 54>C{<j)). 

After this detour, we now turn our attention back to 
the case of interest, the implementation of a coordinate- 
dependent potential. Our discussion is based on the 
precedent one though now generalizing it to the case of 
a coordinate-dependent potential V{x'^,(t>). 

A naive first approach would be to replace occurrences 
of V{(l)) in ^ by y(a;'^,0). However, this will bring an 
unfortunate consequence, namely that one can now no 
longer express the divergence of in the form given by 
cqn ([5]), where Vf,0 appears as a common factor. Instead 
one has 



o = Var% = (Vb0)( v,v'^0-i^F(x^0) 



1 d 



(8) 



The crucial difference with eqn. ([5]) is that several (in- 
dependent) equations must be satisfied by the real scalar 
field (p. As a result, the system of equations will be gener- 
ically inconsistent. 

To resolve this problem we start by: (i) adopting a 
different ansatz for T^P\b (equation ([9]) below), and (ii) 
imposing symmetry conditions on the scalar field. 

First, consider setting T'-P\b, instead of being given by 
equation ([3]) , to be the product of a function of and a 
coordinate dependent tensor. 



(9) 



where the function / is independent of x'^ and the tensor 
Hab is independent of (/). Now, find a suitable Hab such 
that V aT°'b takes the form of equation (O, this will in- 
duce conditions on Hab- Under this choice the divergence 
of the stress-energy tensor results 

df 

(10) 

Now, we look for conditions that would allow us to ex- 
press the r.h.s. of equation pop in such a way that Vb0 
appears as a common factor. Since H'^b is independent 
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of (j), Vb4> cannot appear in the last term of (jlOp: Then, 
that term must be zero, resulting in the first condition 
on 



Va-ff% = 



(11) 



We now consider the second term in the r.h.s.; the con- 
dition 



(12) 



for some scalar h(x'^), ensures that that term has Vb0 
as a common factor. Equation ()12p is satisfied for any 
scalar field 6 if 



(13) 



However, this condition, together with equation im- 
plies that h{x'^) is a constant. This means that T'^P^b is 
of the form ([3]) (with V independent of x'^). Thus, for an 
arbitrary scalar field, and without any further structure 
in the spacetime, space-dependent potentials can not be 
considered. However, by imposing further conditions on 
the scalar field (/>, H"'^ can indeed be chosen with further 
structure than that of equation while still satisfying 
equation (fT2)l . To this end, we consider 3J| the tensor 
_ff°f, of the form 



H\ h{x)S\ + A\ 
Replacing ([H]) into we find 

{Vacp)A\ - . 



(14) 



(15) 



The simplest case is the one with A°-b — ^ for which 
H°'b is given by ([13]). More general cases arise when is 
independent on one of the coordinates, lets say d^scp = 
Vac/) — 0. Here one can adopt A'^^ arbitrarily and set all 
other components to zero, thus satisfying equation p^ . 
In this particular case, takes the form 



/i 

/i 

Q {) h {) 

6 



(16) 



for some functions h{x'^) and b{x'^). 

Similarly, when does not depend on two of the coor- 
dinates, lets say d^^cj) = 0, dx3(t> — 0, one can choose 



/i 

/i 

6 

c 



(17) 



Analogous results are obtained when some of its deriva- 
tives are linearly related. For example, if dx!i(f> = cd^^cl), 
one can adopt A^^ arbitrarily and set A^^ — — 
keeping all other components zero. With this choice, 
equation will be satisfied and H'^b will then be given 
in terms of two functions h{x'^) and b{x'^) in a slightly 



different way as is (|16p . 

Summarizing, we have seen that a coordinate depen- 
dent potential can be implemented if the following condi- 
tions are satisfied: (i) its derivatives are linearly depen- 
dent (this includes the possibility of one or more of them 
being zero) . (ii) The "potential" part of the stress-energy 
tensor is given by with satisfying Va-ff°6 = 
and being expressible in the form (fT6|) . (fT7|) . or similar 
expressions depending on how condition (i) is fulfilled. 

In the next section we will consider in detail the case 
of spherical symmetry. 



A. Spherical Symmetry 

We will now concentrate on the case of spherical sym- 
metry. The line element can be written in the form 



-N'^dt^ + grr{dr + (idtf + gndVl^ 



(18) 



where N , grr, /3, and gn are functions of t and r. We 
assume that we can adopt coordinates so that 



difCj) — 0. Then, H°-b is given by pT|) . with the additional 
condition that h = c due to the spherical symmetry. H°'b 
is then 



/i 

/i 

6 

6 



(19) 



with h and 6 functions of t and r. 

The evaluation of V aH'^b gives rise to non-trivial equa- 
tions only on the t and r components, 



dan , , dh 
Jii(/.-6) + 2,.-.0, 

'^ih-b) + 2gJJ^ = 0. 
dr dr 



(20) 
(21) 



In order to obtain a family of solutions to these equa- 
tions we will demand that h depends on the coordinates 
only through go,: h{t,r) — h{gfi{t,r)). With this condi- 
tion, we have that 



dh 

dx^ 



dh dgn 
dgn dx"- 



(22) 



for a;* — {t, r). Substituting this into either equation ((20|) 
or (I21|) . we obtain an expression for 6 in terms of /i. 



9n 



dh 
dgn 



(23) 



We have just seen that, if h depends on the coordinates 
only through go, and 6 is given in terms of h by (|23p . the 
prescription p9|) for the tensor H°'b allows us to express 
^ aT'^b with Vb4> as a common factor. More explicitly: 



df ^ 



(24) 
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Notice that, if one wanted to calculate VaT^f, without 
setting dgif) = d^(f> = at the onset, one would obtain 
(p4|) . but with h{gn) replaced by b{gn) for the angular 
components V a.T'^e and V aT°'Lp- However, because those 
terms are actually multiplied by zero, equation (|24p is 
true for all four components. 

Setting the r.h.s of ([M)) to zero we obtain the equation 
of motion for 0, 

V,VV+|^/i(5o) = 0, (25) 

where we remind the reader that / is an arbitrary func- 
tion of (j), and h is an arbitrary function of g^. 

Throughout the rest of this work we will choose these 
functions as 

/(</') = -\<l^\ (26) 
h{gn) = + Vign). (27) 

We do that, so that the equation of motion for the scalar 
field becomes 

(VaV^ - m2 - V{gn)) = 0, (28) 

where we interpret the function V{gfi{t,r)) as a 
(coordinate-dependent) potential. The parameter m is 
set to zero in our simulations. The function gQ,{t,r) is 
just the square of the areal radius, R{t, r). Then, we can 
write in the form 

(VaV" - - V{R)^ = 0, (29) 

where V is an arbitrary function of the areal radius. 

In appendix we summarize the results obtained in 
the case of axial symmetry. 



III. THE EQUATIONS 

In this work we solve the non-vacuum Einstein equa- 
tions for a dynamic spherically symmetric space time, 
coupled to a real scalar field. The scalar field satisfies a 
Klein-Gordon-like equation with the addition of a poten- 
tial, as explained in section Hi Al 

The equations are decomposed using a Cauchy formu- 
lation, in which the space-time is foliated by space-like 
surfaces. The particular formulation used is the Einstein- 
Christoffel hyperbolic formulation jl7j, where the equa- 
tions are decomposed into a system of first order hyper- 
bolic "evolution equations," plus a system of (first order) 
"constraint equations." These equations can be solved by 
giving initial data that satisfy the constraint equations on 
a given surface of the foliation, and then integrating the 
evolution equations in time. The constraint equations at 
later times are then automatically satisfied 16] in the 
domain of dependence of that surface. 



The equations solved are the Einstein-Klein-Gordon 
equations, with the addition of a potential. 

Gab = SirTab , (30) 
(VaV'^-y)</. = , (31) 

where the stress-energy tensor. Tab, and the potential, V, 
are given according to section III Al as well as the condi- 
tion that is independent of {6, 0). In equation pip we 
have set m = 0, but this parameter can be incorporated 
in the definition of V. 

We consider the line element and extrinsic curvature 
of a space time in spherical symmetry in the form 

ds^ = -N^dt^ + grr{dr (idtf + r^gTdnl^2) 
K.jdx'dx^ = Krrdr^ +r^KTdfl^, (33) 

where /3 is the (r component of the) shift vector, and 
N is the lapse function. In the Einstein-Christoffel 
formulation, the shift and "densitized lapse" function, 
a = N/ y/g, are arbitrarily specified and kept fixed dur- 
ing the evolution. We denote by g the determinant of the 
three-metric. 

In spherical symmetry, this system reduces to nine first 
order evolution equations, and four first order constraint 
equations, the later containing only spatial derivatives. 

The variables evolved are: the metric components, grr 
and gx', the scalar field, 0; and other variables used to 
convert the equations from second to first order. They 
are: the extrinsic curvature components, Krr and Kt 
(defined in eqn. ([551) ): variables {^,11} constructed with 
first-derivatives of 0, 

* = a,0, (34) 

n = ^(/3 9.0-9t0); (35) 

and the variables {frm frx} containing first spatial 
derivatives of the metric, 

drgrr , '^9rrfrT 
Jrrr = h , (36) 

2 gr 

/.T = %^ + ^. (37) 
2 r 

The complete expressions of these equations are shown 
in detail in appendix IB| Their derivation, and the nota- 
tion used, is based on |l8| and p^ . with the addition of 
terms containing the potential. 

IV. NUMERICAL IMPLEMENTATION 

A. Initial Data 

Consistent initial data must satisfy equations (|B12[) - 
(IB16p . These equations determine some variables in 
terms of others judiciously chosen. In this work, we 
exploit this freedom to describe a black hole centered 
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at r = by specifying {V, (p, gm Krr} from the known 
Schwarzschild solution and solving for gj- and Kt- 

Before describing the details of our implementation, we 
discuss how the potential and scalar field are chosen. We 
adopt a potential V with two free parameters {A, tq} to 
regulate the depth and location of the "well" where the 
scalar field is to be confined (see figure [1]). A simple 
expression for V suffices for this task, and we adopt 



V{R) = A ( 



1 



-{R-raY 



(38) 



with the areal radius R given by i? T-Jgr- The pa- 
rameters in this expression were set \,q A — and 
ro = 6M, where M is the initial mass of the black hole. 
Notice that during the evolution R = R{t,r), thus, in 
these coordinates, the shape (and position) of the poten- 
tial well can change in time. We will return to this point 
later. 

The scalar field (p is defined following cither one of 
two different strategies. One is designed to conform to 
time-harmonic situations in weakly-gravitating cases and 
the other simply prescribing a sufficiently smooth pro- 
file. The latter choice allows us to investigate the space- 
time's response to fields not designed to conform to a 
time-harmonic dependence. 



1. Time-harmonic scalar field 

To prescribe a scalar field which will give rise to a 
spacetime with harmonic time-dependence, we begin by 
considering the limiting case when the scalar field's am- 
plitude is negligible; there the metric should be described 
by the Schwarzschild's solution. Now, considering the 
scalar field as existing over this fixed background space- 
time, a Schrodinger-like eigenvalue equation can be ob- 
tained to determine time-harmonic states as discussed 
below. 

The Schwarzschild metric in Eddington-Finkclstein co- 
ordinates is: 



r / 



2M 



dr' 



AM 



dtdr + r'^dn'^. 



(39) 



We use this metric to evaluate the equation of motion for 
(j>, equation ([31]) . To solve this PDE we use the following 
ansatz that yields separation of variables [35], 



4){t, r) — u{r) cos oj 



t - 2M In 



2M 



M 



(40) 



The equation for u{r) results 



C u{r) 



r / 



V{r) 



u{r), (41) 



where the second order operator C is given by 
£ = - I 1 



2My 5^ 
r I dr^ 



1 



M 

r 



1 



2M 

r 



d_ 

dr 



(42) 



Equation (|4ip is integrated to obtain u{r). Then, from 
its definition, equation P0|) . (/)(t, r) is calculated. Finally, 
from (pit, r) we obtain Ii{t, r) and $(t, r) evaluating these 
functions at t = and adopting them as initial data. 

Equation (|¥T|) can be straightforwardly integrated to 
obtain both the eigenvalue and eigenfunction through a 
standard shooting algorithm. To this end, we transform 
the second-order equation to a system of two first order 
equations for u(r) and u'{r) = du/dr augmented with a 
third equation ((w^)' = 0) to simplify the implementation 
(see (2^] for the details). 

The system of equations is then integrated outwards 
from = 4M on one hand, and also inwards from 
rji = 8M. The obtained solutions are matched at an 
intermediate point, in our case at vq (the center of the 
potential well), with the conditions that both the solu- 
tions and derivatives are continuous. The initial guesses 
for the boundary conditions are then varied until a satis- 
factory match is obtained. The code used to implement 
the shooting algorithm is the one described in [l^] , except 
that the ODE integrator is replaced for LSODE (Liver- 
more Solver for ODEs) |21j]. The boundary conditions, 
consistent with the physical scenario in mind are deter- 
mined as follows. 

We have a system of three first order ODEs, thus three 
boundary conditions need be specified. Natural condi- 
tions for our purposes result from requiring the fields fall 
sufficiently rapid at the boundaries. We thus impose a re- 
lationship between u and its derivative at each boundary, 
of the form u' = ku. The coefficient k at each boundary 
can be found through a WKB-type approach. To do so, 
we first consider the variable change u{r) = F{r)u{r) and 
fix F{r) = [r{r — 2M)]~2 so as to remove the first order 
derivative in equation (|4ip . The resulting equation is 



/(r)u"(r) + Veff{r)u{r) = uj^u{r) 



with /(r) and Vcs{r), 

fir) = (l 
V,tf{r) = ( 1 



2M 

r 
2M 

r 



M2 

V{r) —. 



(43) 

(44) 
(45) 



and we interpret VcS as an effective potential (which is 
shown in figure [T]). Next, we freeze the coefficients /(?■) 
and Vcs on a small neighborhood of each boundary point 
and consider solutions of the form exp(±/cr), with fc^ = 
(V^ff — o;^)//. The following conditions at the boundaries 
are then determined by 

u{r) (X e+'^i'^ 



at 



u{r) (X e 



-k-2r 



at r = tr, 



(46) 
(47) 
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FIG. 1: Potential and effective potential, as defined in (|38p 
and (|45[) . respectively. As mentioned in the text, the poten- 
tials are, in general, functions oi R = r^fg^. The potentials 
showed in this figure are those used to find the time harmonic 
states «(r), where the Schwarzschild metric is used, hence 
R = r. 



where 



'l4ff(r-L)-w2 



(48) 
(49) 



As illustrated later, these conditions indeed ensure the 
solutions decay rapidly outside of the potential well (for 
a bounded range of values of w^). Notice that since the 
equations are homogeneous there remains a freedom on 
the amplitude of the fields at the boundaries. We fix 
this freedom by setting m = I at and adopting as the 
varying parameter for the shooting method the value of 
w at rfl. 

Once obtained (^(r) in [ri,ri^] using equation (HD|) . we 
set ^{t) = outside this region. For the amplitude of 

used in this work in the case of time-harmonic initial 
data, the values of ^ and its derivative at and rj^ are 
small enough to ensure that this matching is sufficiently 
smooth, as is corroborated when evolving these initial 
data. 



2. Smooth profile 

The other approach employed in this work is to adopt 
a simple expression for the scalar field. In particular we 
adopt a "pulse" of compact support of the form 



(/.(r) = 



c{r — ri)*(r — r2)^ fi < r < r2 







elsewhere 



(50) 



where the values ri and r2 control the width of the pulse 
and were chosen so that it is centered with the potential 



(at r = 6M): ri = 5M, r2 = 7M. After specifying ri 
and r2, the coefficient c is chosen so that the scalar field 
has a given mass. This initial data is used to compare 
with the previous approach in regimes where the fixed- 
background approximation is justified and to study the 
spacetime's behavior in non-linear cases. 



Remaining data 

Having specified both the potential and the scalar field, 
consistent initial data is determined by integrating the 
constraint equations in the following manner. First, the 
functions grr, Krr, Oi, and /3 are set equal to those read-off 
from the Schwarzschild solution in Eddington-Finkelstein 
coordinates. Adopting these coordinates gives the free- 
dom to place the inner boundary inside the black hole. 
Wc found it convenient to rewrite the constraint equa- 
tions in the form: 



drOT = dx: 

drdr = fi{gT,dT,KT;Fi), 
drKr = f2{gT,dT,KT;Fi), 



(51) 
(52) 
(53) 



where Fi represents all the functions that are specified 
a priori (including 0). These equations are integrated 
outwards from the inner boundary using the step adap- 
tive integrator LSODE, using as boundary data {gr, 
dx, and Kt at r ~ rmin) the values read-off from the 
Schwarzschild solution. 



B. Evolution 

We discretize the equations with a scheme formulated 
to take advantage of numerical techniques which guaran- 
tee stability of generic linear first order hyperbolic sys- 
tems. In this work we adopt: (i) second order accu- 
racy by implementing second-order derivative operators 
satisfying summation by parts [H, HE 113, HE Ell ; (h) 
a third-order Runge-Kutta operator for the time inte- 
gration through the method of lines ^7}; (iii) a Kreiss- 
Oliger [11] style dissipative algorithm to control the high 
frequency modes of the solution [H,[l^,[33| and (iv) max- 
imally dissipative boundary conditions setting all incom- 
ing modes to zero [1^ Hlj . 

We employ a uniform grid to cover the region r S 
[I'min, ''max] with N equi-spaced points. The grid-spacing 
between points is Ar — (rmax — '"min)/(-/V — f ). The time 
step At is defined in terms of Ar as At = cfl Ar and 
cfl = 0.25 is chosen so that the CFL condition [32] is 
satisfied. In what follows, sub-indices denote particular 
points of a slice, and super-indices distinguish each slice. 

The inner boundary, r = r-min, is set inside the black 
hole initially, and monitored during the evolution to en- 
sure that it remains inside and constitutes and outfiow 
boundary of the computational domain. Then, there is 
no need to prescribe boundary conditions there. At the 
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outer boundary, r = Tmax maximally dissipative bound- 
ary conditions are adopted. In our present case we take 
the simplest form of these conditions and set the incom- 
ing modes to zero. The characteristic structure for the 
system of equations is detailed in appendix [Cl 

The code have been tested to ensure that the numeri- 
cal solutions obtained converge to the corresponding so- 
lutions of the Einstein equations. In appendix [P] we show 
the convergence test for the Hamiltonian constraint. 

V. ANALYSIS AND RESULTS 



space-time. By varying the initial guess for the frequency 
in the shooting integration we obtain different modes. 
We show the first modes in figures [2] and [3l However, for 
this work we used only the first mode which will be re- 
ferred to as "the timc-harmonic state" , unless otherwise 
specified. These modes have been re-scaled so that they 
can be normalized (in analogy with quantum mechanics) 
so that J r'^\u{r)\'^dr = 1. There is no physical justifica- 
tion for choosing that particular normalization, but it is 
helpful when comparing different eigenstates, which oth- 
erwise would have greatly different amplitudes. 



In the simulations performed in this work we set the 
initial mass of the black hole to M = 1 (in geometrized 
units). The domain of integration was chosen so that the 
region of interest is unaffected by the conditions adopted 
at the right boundary. This corresponds to ry^in = 
and Tmax = 221Af . The maximum resolution used was 
Ar = O.OIM (22000 grid points). 

In the two approaches we use to obtain initial data, we 
have the freedom of adjusting the amplitude of the scalar 
field, which in turn determines its mass. We set initial 
data where the mass of the scalar field is irist = O.OIM in 
the time-harmonic case, while for the non-time-harmonic 
cases we set mgf equal to O.OIM, kO.IM, (M being the 
initial mass of the black hole and k = 1...5). To calculate 
the mass we use the Misner-Sharp formula 



MMs(r) 



1 



9T 



f2 



(54) 



which measures the total mass inside a spherical surface 
labeled by coordinate r. In our initial data the mass of 
the black hole, M, is preset, so we can calculate TOsf by 
subtracting M from the total mass of the space-time, 



(55) 



where R labels a sphere containing the scalar field, which 
is localized initially. (See figure [8]). 

During the evolution we employ this formula, replacing 
M for Mms at the horizon |36l|. 

In our analysis we also evaluate the Kretschmann in- 
variant I = Rabcd.R"''^"^ , where Rabcd is the Riemann ten- 
sor. This quantity provides a gauge-invariant answer that 
can be compared with its value in known spacetimes. For 
a Schwarzschild space-time, / is given by 



'Sch 



48(Mms)' 

i?6 



(56) 



where, in Schwarzschild coordinates, Mms — M and R = 
r. We evaluate the quotient I/Isch using with R = 
r^/gr and Mms defined in ((54)) . 

A. Initial Data 



As explained in section IIV A 1[ we first find time- 
harmonic states for the scalar field on a Schwarzschild 



first mode 
second mode 
third mode 




I 

6 
r/M 



FIG. 2: First time-liarmonic states of u{r). 
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-2e-03 
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— first mode 

- - second mode 

third mode 




FIG. 3: Scalar field at t = obtained from the first time- 
harmonic states of u(r), using equation (|40p . 



The other approach used to define the initial data cor- 
responds to the "pulse" described in section [TV A 21 In 
the linear regime we employ both types of initial data, 
with a scalar field's initial mass nisf — O.OlAf . In the 
non-linear regime we adopt only the non-time-harmonic 
initial data with masses mgf ranging from O.IM to 0.5M. 
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B. Evolution 

We study the evolution of the prescribed data. We be- 
gin by considering first the hnear regime, adopting scalar 
field configurations with initial mass of 1% of that of 
the black hole. After confirming that the time-harmonic 
configuration behaves as expected, we confirm that the 
"pulse" configuration evolves towards a time-harmonic 
regime. Then, we study cases in the non-linear regime, 
with initial scalar field masses ranging from 10% to 50% 
of that of the black hole. In all cases we evolve until 
t = 200M. 



about 10% of the field's initial mass falls into the black 
hole. 

The amount of mass that falls into the black hole is 
calculated by subtracting the Misner-Sharp mass at the 
horizon, minus the initial mass of the black hole. In the 
case of time-harmonic initial data this number is (1 ± 
3) X 10~^Af, while for that of non-time-harmonic initial 
data it is (10 ± 3) x IQ-^M (see table Hand figure [HI). 
These values are calculated using the highest resolution 
(Ar = 1/lOOM), and the errors as the difference of these 
values with those of a lower resolution (Ar — 1/50M). 



Linear case 

The time-harmonic initial data constructed essentially 
remains unchanged through the evolution while the non- 
time-harmonic data evolves towards a time-harmonic 
state. Figures|4]and[5]illustrate 0(r) at different times for 
the maximum resolution employed (Ar — M/lQi)). Fig- 
ure 2] corresponds to the time- harmonic initial data, and 
figure [S] to non-time-harmonic initial data. In both cases 
we sampled along two different periods aX t k, 80Af ; and 
then at t K, 160M. The corresponding pairs, are then 
plot together illustrating how after 22 periods apart the 
solutions are essentially the same. This is further illus- 
trated in figure [6] where we show the difference between 
each of these pairs for three different resolutions. 

Finally, figure [7] displays the absolute value of the 
Fourier transform in time of ^(j) dr, denoted The 
scalar field is first integrated in space, then a discrete 
Fourier transform in t is calculated, where t ranges from 
to 200M in the case of time-harmonic initial data, and 
from to = 60M to 200M in the non-time-harmonic case. 
In the plot we also indicate the frequencies (/„ = w„/27r) 
obtained from the shooting integration when calculating 
the time-harmonic states. The time is chosen after the 
initial transient behavior, indicated by a time-harmonic 
behavior observed in (j). The initially non-time-harmonic 
scalar field relaxes to a superposition of the first three 
time-harmonic modes, the first one being the dominant 
one. We point out here that for this configuration, the 
shooting method gives raise to three possible modes. It 
is thus no surprising that the evolution gives rise to a so- 
lution described by these modes. Deeper potentials give 
rise to more modes. 

Figures [5] and ini show the Misner-Sharp mass function 
(equation (HH)) for both types of initial data. The con- 
tinuous line shows the initial value (Mms at t = 0). The 
discontinuous lines show A^ms at t = 200M for three dif- 
ferent resolutions. In both cases the asymptotic value 
of the mass stays constant, indicating no scalar field en- 
ergy is radiated away. An inspection of the mass behav- 
ior at smaller radii for the solution obtained with time- 
harmonic initial data reveals that this converges to essen- 
tially the initial value, thus a negligible amount of mass 
falls into the black hole. For the non-time-harmonic case 



Non-linear case 

We turn now to the non-linear cases investigated. 
These correspond to initial mass configurations where the 
scalar field has a mass of at least 10% of that of the black 
hole. In this regime we solely adopt the "pulse" prescrip- 
tion defined in equation ([50]) for the scalar field since 
the time-harmonic data is obtained under an assumption 
which is no longer valid. 

As we have done for the linear case, we also compare 
profiles at different times for simulations with higher ini- 
tial TTisf. Figures [TUl and [TT] correspond to initial masses 
of the scalar field of nisi = O.lOAf and rrist — 0.50M, 
respectively. The time it takes to reach a state described 
by a harmonic time dependence is longer than in the lin- 
ear regime, especially for the higher initial rrisf — O.SOAf . 
For that reason, the first samplings (labeld ti in the fig- 
ures) occur later than in the linear case, and the interval 
between the profiles compared, ^2 ~ ii, is ten periods, as 
opposed to 22 in the hnear cases. 

The absolute value of the Fourier transform oi Jcj) dr, 
\F[(I)]\, is shown in figure [T^ for the two different initial 
masses of 0. Again, we compute the transformation after 
the initial transient behavior has passed and the scalar 
filed has already reached a quiescent state. As a useful 
indicator, we also show the frequencies corresponding to 
time-harmonic states. Now, while the observed modes 
do not coincide exactly with those obtained at the linear 
approximation, they are close to them. 

In figure [T^ we show the Misner-Sharp mass at t = 0; 
and a,t t = 200Af for three different resolutions. Fig- 
ures 13(a) and |13(b)] correspond to initial masses of the 
scalar field of mgf = O.lOAf and mgf — O.SOAf, respec- 
tively. In all these cases about 10% of the scalar field's 
mass falls into the black hole, while nothing escapes out- 
wards. Additionally, for the case with grater mass, the 
scalar filed spreads slightly outwards before reaching a 
quiescent state. Although we only show figures corre- 
sponding to two different initial values of mgf, we have 
simulated the system for other values of this parameter 
nisf — k10~^M (k — 1...5). In all these cases essentially 
no scalar field energy is radiated away, while a small por- 
tion falls into the black hole. The measured values are 
shown in table U and figure [TH 

If, after some transient time, the scalar field is finally 
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r/M 



(b) 



FIG. 4: The scalar field at different times is compared to check if the evolution remains described by a time-harmonic depen- 
dence. Case with time-harmonic initial d ata. I nitial mass of the scalar field msf — O.OIA/. Figure [4(a) | shows the scalar field 
when it reaches a maximum, while figure 4(b) shows it at about a quarter of a period later. In both cases, the profile shown 



in continuous line is separated 22 periods from the one in dashed line. 




FIG. 5: Here we show the same comparison of profiles as in figure|4l this time for the case with non-time-harmonic initial data. 
The separation between the profiles compared is also 22 periods. The initial mass of the scalar filed is rrisi = Q.OlAf . 



TABLE I: Mass that falls into the black hole for different 
initial masses of the scalar field. Calculated as the Misner- 
Sharp mass at the horizon ai t — 200 minus the initial mass 
of the black hole. See figure 1141 

Initial m^f [M] (MMs(nO - M) [lO'^M] 

0.01 0.10 ±0.03 

0.10 1.0 ±0.3 

0.20 2.9 ± 0.7 

0.30 3 ± 1 

0.40 5 ± 1 



confined within a compact region, lets say [ra,rb], the 
space-time should be that of Schwarzschild for r > r^, 
with a Schwarzschild mass equal to the total mass inside 
the sphere r = r?,. This can be checked by evaluating the 
Kretschmann invariant. In figure[Tn]we show the quotient 
///gch (see the paragraph containing equation (|56p ) at 
t sa 140M for the case with initial mgt = 0.5M. This 
quotient converges to one for r > , and also for r < . 
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(a)|0(ti) — </'(t2)l for the case with time-harmonic initial data (see {h)\<f>{ti) — 4'it2)\ for the case with non-time-harmonic initial data 
figure [T( a) ^ (see figure [5'( a) | 



FIG. 6: Absolute value of the difference between the scalar field at different times: \<j){ti) ^ 0(^2)|, where t2 — ti — 22 periods. 
Figure 6(a) shows the difference between the profiles shown in figure 4(a) while figure 6(b) shows the difference between those 
in figure |5(a)] In each case, we show these differences for three resolutions. 
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FIG. 7; Absolute value of the discrete Fourier transform in 
time oi j(f> dr. The continuous line corresponds to the time- 
harmonic initial data, while the dashed line corresponds to the 
non-time-harmonic initial data. In the later case, the scalar 
field relaxes to a superposition of the first time-harmonic 
modes, whose frequencies are shown in the figure (labeled 

/n). 
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FIG. 9: As in fig. |8l we show the mass function, this time for 
the non-time-harmonic case. The initial mass of the scalar 
field is nisi = O.OlAf. This time about 10% of it falls into the 
black hole. 
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FIG. 8: Mass function at t = 0; and &t t = 200M for three 
resolutions. Stationary initial data. Initial rrisf = O.OIM. The 
continuous line shows the mass function at t = 0, while the 
discontinuous lines show, for difi^erent resolutions, the mass 
function a± t = 200M. In this case the escape of mass into 
the black hole is negligible (Arrisf = (1 ± 3) x lO'^M). 
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FIG. 10: The scalar field at diflterent times is compared to check if the solution obeys a harmonic time dependence. Case with 
non-time-harmonic initial data. Initial mass of the scalar field nisi ~ O.lOAf . Figure 10(a) shows the scalar field when it reaches 
a maximum, while figure [lO(b)| shows it at about a quarter of a period later. In both cases, the profile shown in continuous 
line is separated 10 periods from the one in dashed line. 
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FIG. 11: This figure shows the same comparisons as figure [10] but for an initial mass of the scalar field of rrisf = 0.50M. The 
separation between the profiles compared is also 10 periods. 
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FIG. 12: Absolute value of the discrete Fourier transform in t of the space integral J <j>{r, i) dr. The marks labeled /n denote 
the frequencies of the first modes obtained from the shooting. The three peaks, which indicate the dominant frequencies in the 
solution, lie at slightly lower frequencies than those of the time-harmonic states in the linear case. This behavior is consistent 
with the frequency shift due to the black hole growing in size. However, the growth alone does not fully account for the observed 
shift, though this is expected as non-trivial contribution due to non-linearities also play a role. 
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FIG. 13: Mass function sX t = Q\ and at f = 200M for three resolutions. The discontinuous lines show the mass function at 
t=200M for three resolutions. In each of these cases, about 10% of the initial mass of the scalar field falls into the black hole, 
while nothing escapes to infinity. 
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FIG. 14: Mass that falls into the black hole for different initial 
masses of the scalar field. Calculated as the Misner-Sharp 
mass at the horizon at t = 200 minus the initial mass of the 
black hole. See table U 
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FIG. 15: Kretschmann invariant quotient for three resolutions 
at ti — 139.983Af. This quotient converges to 1 outside of the 
region where the scalar field is confined. A horizontal line at 
I /Isch = 1 have been drawn as a guide. Initial rrisf — 0.50M. 
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VI. CONCLUSIONS 

We have discussed difficulties encountered when at- 
tempting to confine a scalar field distribution within some 
region. The existence of a symmetry in the spacetime al- 
lows for doing so in a consistent manner. For the specific 
spherically symmetric case, we have given prescriptions 
for implementing a scalar field with a potential depend- 
ing on the areal radius R. 

We have illustrated the viability of this approach by 
confining a scalar field distribution around a black hole. 
For our particular choice of potential and initial scalar 
field, the scalar field becomes totally confined after some 
transient time, which depends on the initial mass. During 
the transient, part of the scalar field accretes into the 
black hole, while nothing escapes to infinity. By adjusting 
the depth of the potential, the amount of energy that falls 
in can be controlled. 

The approach can be exploited, an extended, to mimic 
situations of interest. These can range from physical 
studies of particular systems, to serve as a testing model 
for infrastructure development aimed to simulate more 
complex systems. 
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APPENDIX A: COORDINATE-DEPENDENT 
POTENTIAL IN AXIAL SYMMETRY 

Following fsS^, we write the (general) axi-symmetric 
line element in the form 

-e^^idtf + e^'^' {dip - qidx^ - q2dx^ - udtf 



ds' 



(Al) 



where all the functions appearing here are functions of 
= t, a;^, and x^, but independent of x^ = if. We 
assume that the scalar field (f> is independent of ip, and 
hence use 77"^ as given in (fT6|) . Evaluating V aH"-^ (and 
assuming that h and b are independent of ip) we find 
that the (/s-component is of the form [h — b) times an 
expression depending on the metric functions and their 
derivatives. We assume that the expression multiplying 
(h — b) is not zero, because at this moment we want 
to consider the case of no other symmetry other than 
the axial symmetry. Then, setting this component to 



zero, we have the condition b = h, which, as we have seen 
earlier, implies that h is a constant. This means that the 
potential will be independent of the coordinates. 

Consider now the special case of axial symmetry with- 
out rotation. We can write the line element in the form 

ds^ = g^p^dp^ -\- gijdx^dx-' , (A2) 

Evaluating \I aH'^bi the (^-component this time results 
identically zero, and, setting the other components to 
zero, we have: 



dx"^ 



dh 



(A3) 
(A4) 
(A5) 



At this point, one can follow the same procedures as in 
section Hi Al (compare these equations to (|20|) and ((2T|) '). 
For that reason, in this section we will just summarize 
the results. 

Equations (jA3|l - (jA5|l are satisfied if, (i): h depends 
on the coordinates only through an arbitrary function of 



h{t,x^,x'^) = f{g^^{t,x^,x^)), 
and (ii): b is given in terms of h by 

dh 



b = h + 2g^^ 



dg^ 



(A6) 



(A7) 



Given these conditions, one can express 'VaT°'b with Vb</) 
as a common factor, and, setting it to zero, obtain the 
equation of motion for the scalar field. 



df 

VaV> + /i(.9^^) = 0, 



(A8) 



where / and h are arbitrary functions of (p and g^plp^ re- 
spectively. On can, in particular, choose these functions 
as follows, 

1 



(A9) 

/i(ffw) = m^ + U(g^^). (AlO) 

Then, the evolution equation becomes 

(VaV^ - m2 - U{g^^)) = 0, (All) 

where we can interpret C/ as a coordinate-dependent po- 
tential. 

APPENDIX B: THE EQUATIONS 

The equations of motion are 
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9rr = Pg'rr + '^9rrP' " 2agU''gTKrr , (Bl) 

gr = Pg'r - ^dg^grKr + ^ , (B2) 

r 

krr = PK'„ - ag^^^/^grfrrr " ""s'/'st " Ggr^gU^af^T + 4.gTr-^gU^a' - dgrr-^gH^a + 2if„/3' (B3) 

-grgrr^^aKl^ + 2gll^aKrrKT - Sg^r^'^afrTfrrr + 2gTgrr^^0^frrr 
+ 2gTr~^grr^^afrrr " gTgrr^^&frrr + grgW^ Ct Att(T g^r - 2Srr) , 

Kt = PK^ ~ agrg-r'^^frT + 2Pr-'KT + grr-^g^a + agTRrK^g^r^/^ ~ grfrT^'grr'^ - 2af^Tgrr'^ (B4) 

+agll''gT MTgr - 2St) , (B5) 

frrr = Pfrrr " agU'^grKr " 'igrl^a'Kr + 12g:^^ g^l^ frT - igll^aKTfrrr - grgrr^^aRrrfrrr (B6) 

-lOgll^aKrrfrT + ifrrrP' + grrP" " ^'gi/'sT^rr + 2r-^ grgU^ aKrr + Sr^^gU^aKr + ^ag^^ gr iirJr , 

frT = PfrT " &gU^9TK^ + /3' frT " a'gU^grKr + 2g'^^aKTfrT ~ ag-r '^Krf rrrgr + 2r-^PUT , (B7) 

$ = agll'^gTli' - g-^^^agrllfrrr + 2agli^IlfrT + 2r-^agll'^gT^ - a'gH^grTl + , (B8) 

il = pW - g-^l'^agr^' + g-^''^agTl^Krr + 2agll'^I{KT - 4g-//2a$/^T + 2r-^g;^''^agT^ - g'^/'^gr^a' (B9) 

+gU''gTaV<j} , (BIO) 

4> = pel)' -grgU^an, (bu) 



where a = ar^ sin 6* = N/^/g^gT] dots denote deriva- (^^^ ^ - 2/rT , (B15) 

tive with respect to t and primes denote derivatives with 



r 



respect to r; and the "source terms" are defined in equa- 
tions ((BT7 l) - (|B2T|) . C'm = * - 0' , (B16) 

The constraint equations are where the "source" terms are difined as 



grrST 2r^gT grrgr 



47rr = -2V(j)^ + 



5rr 



47rp , (B12) 



r^ffTT;— , (B18) 
ogn 

if;; 2iv:T 

Cr = — + ^ + AnJr = $n , (B19) 

gr rgr Qt 

+AlTJr , (B13) 

o r An{Tgrr-2Srr) = -^.^^ 2^2 ^ 

^ / , ogrrJrT r, r /-oi a\ 

Crrr = ffrr H 2/„r , (B14) ^ ^ dV 

9T -^^g^yg^ , (B20) 

ogn 
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47r (TgT - 2S't) = -grVif)^ 
where — r^gr (see section HTXI . 



(B21) 



APPENDIX C: CHARACTERISTIC STRUCTURE 



zero with order two. In figure [16] we show the Hamilto- 
nian constraint, in the case of the strongest scalar filed 
studied, that with initial mgf = O.bM. 

The evaluation of the constraints is a particularly im- 
portant test in this work, to ensure that the implementa- 
tion of a coordinate dependent potential is not breaking 
the covariance of the theory. 



The characteristic modes and eigenvalues obtained at 
a surface r = const are given by 
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APPENDIX D: CODE TESTS 

The standard code tests have been performed, show- 
ing that all the constraints and residuals converges to 



FIG. 16: L2 norm of the Hamiltonian constraint (equation 
(|B12I) ') for three different resolutions. The measured conver- 
gence results of order two as expected. 
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